
    options(error=stop)

    
    
    getlogop <- function(p0,p1){
        log(p0) + log(p1) - log(1-p0) - log(1-p1)
    }
    
    getlogor <- function(p0,p1){
        -log(p0) + log(p1) + log(1-p0) - log(1-p1)
    }
    
    getlogrr <- function(p0,p1){
        log(p1) - log(p0)
    }
    
    #### Function for checking if two things are equal 
    #### within numerical precision
    same <- function(x,y){isTRUE(all.equal(x,y))}
    
    
    ### Function for finding probabilities from logrr and logop
    ###
    getProbScalarRDiff <- function(atanhrd,logop,weight = 1){
        
        rd = tanh(atanhrd) * weight
        
        if(logop>350){
            if (atanhrd < 0){
                p0 <- 1
                p1 <- p0+rd
            }else{
                p1<-1
                p0 <- p1-rd}   		
        }
        else{ ## not on boundary
            if(same(logop,0)){## logop = 0; solving linear equations 
                p0<-0.5*(1-rd)}
            else{
                p0 <- (-(exp(logop)*(rd-2)-rd)-
                           sqrt((exp(logop)*(rd-2)-rd)^2+4*exp(logop)*(1-rd)*(1-exp(logop))))/(2*(exp(logop)-1))}
            p1 <- p0 + rd}
        return(c(p0,p1))
    }
    
    
    Optimize = function(objective, startpars, thres = 1e-6, max.step = 1000){
        
        Diff = function(x,y) abs(x-y)/(x+thres)
        step = 0;    value.old = diff = thres + 1
        while(diff > thres & value.old > thres & step < max.step){
            step = step + 1
            opt = optim(startpars,objective,control=list(maxit=max.step))
            diff = Diff(opt$value,value.old)
            value.old = opt$value
            startpars = opt$par
            if(MESSAGE & step %% 10 == 0){
                cat("This is the ", step, "th step. The optimum value is ",
                    opt$value," after ",opt$counts[1]," iterations \n",sep="")   
            }
        }
        
        opt = list(par = opt$par, convergence = (step < max.step), 
                   value = opt$value)
        return(opt)
        
    }
    
    
    GetName = function(var) deparse(substitute(var))
    GetVar  = function(varname)  eval(as.symbol(varname))
    
    logit = function(p) log(p/(1-p))
    expit = function(x){
      output = rep(1,length(x))
      index  = !is.na(x) & x <= 100
      e.x    = exp(x[index])
      output[index] = e.x/(1+e.x)
      output[is.na(x)] = NA
      return(output)
    } 
    
    Mean = function(x) mean(x, na.rm=TRUE)
    
    NameVar = function(list,name){
      
      varname1 = strsplit(name,split=".",fixed=TRUE)[[1]][1]
      paste(varname1,"=",list[name][[1]],sep="")
      
    }
    
    NameList = function(x){
      
      sapply(1:length(x), function(i) NameVar(x,names(x)[i]))
      
    }
    
    NameFirst = function(names){
      
      sapply(names,  function(name) strsplit(name,split=".",fixed=TRUE)[[1]][1])
      
    }
    
    ParaValues = function(paras){
      
      ###     Create an array containing all the parameter value combinations
      para.string = paras[[1]]
      for(i in 2:length(paras))   para.string = outer(para.string,paras[[i]],"paste")
      
      ###     Extract numerical values from strings
      Extract = function(name) as.numeric(strsplit(name,split=" ",fixed=TRUE)[[1]])
      para.value = apply(para.string,1:length(dim(para.string)),Extract)          # The first dimension contains the results
      dimnames(para.value)[[1]] = NameFirst(names(paras))
      dim.order  = c(2:length(dim(para.value)),1)
      return(aperm(para.value, dim.order))
      
    }
    
    MyAttach = function(list){  # attach a list to the global environment
      for(i in 1:length(names(list))){
        assign(names(list)[i], list[[i]], envir=.GlobalEnv)
      }
    }
    
    Wrap = function(list,width){
      output = list
      for(i in 1:length(output)){
        output[[i]] = paste(strwrap(output[[i]],width=width), collapse="\n")
      }
      return(output)
    }
    
    
    
    Message = function(paras, para, ptm1){
      para.min = sum(sapply(paras[1:(length(para)-1)], min))
      output   = NULL
      if(sum(para[1:(length(para)-1)]) == para.min){
        if(para["n"] > min(paras$n.cand)){
          output = paste(output, " takes ",round((proc.time() - ptm1)[3],
                                                 3),"s \n",sep = "")  
          ptm1 <<- proc.time()
        }
        output = paste(output, "Sample size ", para["n"], sep = "")
        # params.x = DataTrue(para["n"])
      }
      return(output)
    }
    
    # Print = function()
    
    
    
    
    ##### Rounding functions (start) #####
    
    Round = function(x,num){
      #' Round a number by scientific digits
      #'
      #' Args:
      #'     x: The number to be rounded
      #'     num: the number of scientific digits
      #' Returns:
      #'      The rounded number
      #'
      gsub("\\.$", "", formatC(signif(x,digits=num),digits=num,
                               format="fg", flag="#"))
    }
    
    Round_decimal = function(x,num){
      sprintf(paste("%.",num,"f",sep=""), round(x,num))
    }
    
    Round_min = function(x, num){
      # Round large numbers by scientific digits
      # Round small number by decimal places
      output = Round(x,num)
      output[x<1 & !is.na(x)] = Round_decimal(x[x<1 & !is.na(x)],num)
      return(output)
    }
    
    ##### Rounding functions (end) #####
    
    
    
    
    
    
    
    
  